function [fval, eqm] = eqm_ss_iid_cases_2_to_5(x, par, fun, perfect_reallocation, cali_step)

% -----> Input:
% x:   initial value [B, ALPHA_0]
% par: parameters
% fun: functions to be used
% perfect_reallocation: 1 means the economy has efficient reallocation

% -----> Output
% fval: residuals of equilibrium conditions
% eqm: equilibrium variables
% Some note in the following:
% 1). for model cases 2 to 5, we do not have money
% 2). only one EE for capital and we can express everything in closed forms

B           = x(1);
ALPHA_0     = x(2);


if cali_step == 1
    CHI_k   = par.CHI_k;
    CHI_PI  = par.CHI_PI;
    CHI_q   = x(3);
    XI      = x(4);
else
    CHI_k   = par.CHI_k;
    CHI_PI  = par.CHI_PI;
    CHI_q   = par.CHI_q;
    XI      = par.XI;
end
CHI_PI_k    = CHI_PI;
ALPHA       = par.ALPHA;

a_s         = par.a; 
a_b         = par.a;
%%%% Controling the integration
int_method  = par.int_method;
Rel_Tol     = par.Rel_Tol;
Abs_Tol     = par.Abs_Tol;


Z           = 0; % non-monetary case; the real balance Z = 0;
Z_K         = 0;  



if perfect_reallocation == 1 %%%% always meet the highest-producitivty firm
    
    L         = par.EPS_H_BOUND;
    int_l     = min(max(L / par.THETA - (1 - par.THETA - 2 * CHI_PI) / par.THETA *  fun.h_bound(a_b), fun.l_bound(a_s)), fun.h_bound(a_s)); % effective lower bound
    int_h     = min(max(L / (1 - 2 * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)); % effective upper bound

    q_capital = par.EPS_H_BOUND - par.q_average_prod;
    q_money   = 0; % integration w.r.t. additional payment (could come from capital as payment or money as payment)
    q_liq     = 0;
    q_profit  = 0;
    prob_P    = 0;
    
else %%%% meet arbitrary firms. 

    L         = max(par.EPS_L_BOUND,  - (1 - par.DELTA) * (1 - CHI_q - CHI_k) / (1 - par.tau_k) / B);

    q_capital = integral2(@(eps, eps_t)fun.capital_gain(eps, eps_t, a_s, a_b),...
    fun.l_bound(a_s), min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    @(eps)eps, @(eps)fun.threshold(eps, L, CHI_PI, CHI_PI_k, 1),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

    int_l = min(max(L / par.THETA - (1 - par.THETA - CHI_PI - CHI_PI_k) / par.THETA *  fun.h_bound(a_b), fun.l_bound(a_s)), fun.h_bound(a_s)); % effective lower bound
    int_h = min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)); % effective upper bound

    q_money   = integral2(@(eps_t, eps)fun.money_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    +           integral2(@(eps_t,eps)fun.money_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

    q_liq   = integral2(@(eps_t, eps)fun.liq_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    +         integral2(@(eps_t,eps)fun.liq_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

    %%%% for Holmostrom - Tirol benefit
    q_profit   = integral2(@(eps_t, eps) fun.profit_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    +           integral2(@(eps_t, eps) fun.profit_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

    prob_P   = integral2(@(eps_t,eps)fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
            int_l, min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
            @(eps_t)fun.threshold(eps_t,L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
              + integral2(@(eps_t,eps)fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
            min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)), fun.h_bound(a_s),...
            @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b), 'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);        
         
end


% The residuals of the capital Euler equation
fval(1) = (1 / par.BETA - (1 - par.DELTA)) / (1 - par.tau_k)...
    - B * (par.q_average_prod + ALPHA_0 * ALPHA * (1 - par.THETA) * q_capital...
    + ALPHA_0 * ALPHA * par.THETA * CHI_k * (1 - par.DELTA) * q_money...
    + ALPHA_0 * ALPHA * par.THETA * CHI_PI * q_profit); 


w       = par.ETA / (B / par.A / (1 - par.ETA))^((1 - par.ETA) / par.ETA);

c_star  = (XI / (1 - par.tau_h) / w)^(1 / (-par.SIGMA));

K_bar_K = par.q_average_prod + ALPHA_0 * ALPHA * q_capital + ALPHA_0 * ALPHA * q_liq; 
    
% K       = (c_star + par.G) / (B * K_bar_K / (1 - par.ETA) - par.DELTA);
% 
% %Y       = c_star + par.DELTA * K + par.G;
% 
% Y       = B * K_bar_K  / (1 - par.ETA) * K;

if cali_step ==1
    K      = c_star / (B * K_bar_K / (1 - par.ETA)* (1 - par.G_Y) - par.DELTA);
    Y      = B * K_bar_K / (1 - par.ETA) * K;
    G      = Y * par.G_Y;
else
    G      = par.G;
    K      = (c_star + par.G) / ( B * K_bar_K / (1 - par.ETA) - par.DELTA);
    Y      = B * K_bar_K  / (1 - par.ETA) * K;
end
    
H       =  (Y / (par.A * K_bar_K * K)^(1 - par.ETA))^(1 / par.ETA);

E_a_K   = ALPHA_0 * ALPHA * integral2(@(eps_t, eps)fun.reall(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, 1),...
        fun.l_bound(a_s), min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)), ...
        @(eps_t)eps_t, @(eps_t)min(fun.h_bound(a_b), fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1)) ,'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

E_p_K   = ALPHA_0 * ALPHA *...
        ( integral2(@(eps_t,eps)fun.reall(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, 1),...
        int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
        + integral2(@(eps_t,eps)fun.reall(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, 1),...
        int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b), 'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol));       

prob_A   = 2 * integral2(@(eps_t, eps)fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
        fun.l_bound(a_s), min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
        @(eps_t) eps_t, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), 'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);        
    
average_2nd_price = 2 * integral2(@(eps_t, eps) ((1 - par.tau_k) * B * ((1 - par.THETA) * eps + par.THETA * eps_t)...
        + 1 - par.DELTA) .* fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
        fun.l_bound(a_s), fun.h_bound(a_s), @(eps_t)eps_t, fun.h_bound(a_b));    

threshold_iota = ALPHA * (1 - par.tau_k) * B  * par.THETA * q_money;


%%%% using debt first
debt = integral2(@(eps, eps_t)fun.debt_fun(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_k, 0),...
    fun.l_bound(a_b), fun.h_bound(a_b), @(eps)eps, fun.h_bound(a_s), 'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);


weight     = ALPHA_0 * ALPHA;
moment_1   = (1 - weight) * exp(par.log_eps_mu + par.log_eps_sig^2 / 2) +...
    weight * integral(@(x) x .* fun.g(x, L, par.a), fun.l_bound(par.a), fun.h_bound(par.a));

moment_2   = (1 - weight) * ((exp(par.log_eps_sig^2) - 1) * exp(2 * par.log_eps_mu + par.log_eps_sig^2) + + (exp(par.log_eps_mu + par.log_eps_sig^2 / 2))^2)+ ...
    weight * integral(@(x) x.^2 .* fun.g(x, L, par.a), fun.l_bound(par.a), fun.h_bound(par.a));
prod_std   = sqrt(moment_2 - moment_1^2) / moment_1;


%%%5 endogenous entry
        
surplus_entry = ALPHA * integral2(@(eps_t, eps)fun.surplus_buy(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z/K, 1),...
                fun.l_bound(a_s), fun.h_bound(a_s), @(eps_t)eps_t, fun.h_bound(a_b),...
                'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol) * K / w * XI;

surplus_net = max(0.00000001, surplus_entry - par.iota * Z / w * XI);

if cali_step ==1
    fval(2) = E_p_K - par.P_share_target * (E_p_K + E_a_K);
    fval(3) = (E_a_K + E_p_K) - par.R_ratio_target * (E_a_K + E_p_K + par.DELTA);    
    fval(4) = H / par.H_target - 1;
else
    fval(2) = ALPHA_0 - 0.5 * (1 + erf((log(surplus_net) - par.c_mu) / par.c_sig / sqrt(2)));
end

fval = fval * 1000;

eqm.B       = B;
eqm.ALPHA_0 = ALPHA_0;
eqm.L       = L;
eqm.w       = w;
eqm.c_star  = c_star;
eqm.K       = K;
eqm.K_bar_K = K_bar_K;
eqm.Z       = Z;
eqm.Y       = Y;
eqm.G       = G;
eqm.H       = H;
eqm.E_a_K   = E_a_K;
eqm.E_p_K   = E_p_K;
eqm.q_capital = q_capital;
eqm.q_money   = q_money;
eqm.q_profit  = q_profit;
eqm.q_liq     = q_liq;
eqm.prob_A    = prob_A;
eqm.prob_P    = prob_P;
eqm.threshold_iota = threshold_iota;
eqm.average_2nd_price = average_2nd_price;
eqm.prod_std = prod_std;
eqm.debt = debt;
eqm.surplus_entry = surplus_entry;
eqm.surplus_net = surplus_net;

end